/*Algoritmo de Borwein para calculo de Pi usando threads*/

#include <stdio.h>
#include <gmp.h>
#include <semaphore.h>
#include <pthread.h>
#include <time.h>

void *numerador(void*);
void *denominador(void*);
void *func1(void*);
void *func2(void*);

pthread_t th_num, th_den, th_1, th_2;
sem_t sem_num, sem_den, sem_1, sem_2;

mpf_t a, y, aux1, aux2, pi;
int k = 0;
#define PRECISAO 35000000
int main(int argc,char *argv[]){
	time_t ini;

	sem_init(&sem_num, 0, 0);
	sem_init(&sem_den, 0, 0);
	sem_init(&sem_1, 0, 0);
	sem_init(&sem_2, 0, 0);

	mpf_init2(a, PRECISAO);
	mpf_init2(y, PRECISAO);
	mpf_init2(aux1, PRECISAO);
	mpf_init2(aux2, PRECISAO);
	mpf_init2(pi, PRECISAO);

/*Condicoes iniciais do algoritmo*/
	/*a0 = 6-4sqrt(2)*/
	mpf_sqrt_ui(aux1, 2);
	mpf_mul_ui(aux1, aux1, 4);
	mpf_ui_sub(a, 6, aux1);
     
	/*y0 = sqrt(2)-1*/
	mpf_sqrt_ui(aux1, 2);
	mpf_sub_ui(y, aux1, 1);
	ini = time(0);
	for(k = 0; k < 12; k++){
		pthread_create(&th_num, NULL, numerador, NULL);
		pthread_create(&th_den, NULL, denominador, NULL);
		
		sem_wait(&sem_num);
		sem_wait(&sem_den);

		/*So executa qndo th_num e th_den terminarem*/
		mpf_div(y, aux1, aux2); /*y(k+1)*/
		
		pthread_create(&th_1, NULL, func1, NULL);
		pthread_create(&th_2, NULL, func2, NULL);
		
		/*So executa qndo th_1 e th_2 terminarem*/
		sem_wait(&sem_1);
		sem_wait(&sem_2);

		mpf_sub(a, a, aux2); /*a(k+1)*/
	
		mpf_ui_div(pi, 1, a);
		gmp_printf("PI(%d): %.10000000Ff\n", k, pi);
	}
	
	printf("tempo: %d segundos\n", time(0)-ini);
	mpf_clear(a);
	mpf_clear(aux1);
	mpf_clear(aux2);
	mpf_clear(pi);
	mpf_clear(y);

	sem_destroy(&sem_num);
	sem_destroy(&sem_den);
	sem_destroy(&sem_1);
	sem_destroy(&sem_2);

	return 0;

}

void *numerador(void *arg){
	mpf_t aux3;
	mpf_init2(aux3, PRECISAO);
	mpf_set(aux3, y);/*yk*/
	mpf_pow_ui(aux1, aux3, 4);/*yk^4*/
	mpf_ui_sub(aux1, 1, aux1); /*1-yk^4*/
	mpf_sqrt(aux1, aux1);
	mpf_sqrt(aux1, aux1);/*(1-yk^4)^(1/4)*/
	mpf_ui_sub(aux1, 1, aux1);/*1-(1-yk^4)^(1/4)*/

	mpf_clear(aux3);

	sem_post(&sem_num);
	return NULL;
}

void *denominador(void *arg){
	mpf_set(aux2, y);/*yk*/
	mpf_pow_ui(aux2, aux2, 4);/*yk^4*/
	mpf_ui_sub(aux2, 1, aux2); /*1-yk^4*/
	mpf_sqrt(aux2, aux2);
	mpf_sqrt(aux2, aux2); /*(1-yk^4)^(1/4)*/
	mpf_add_ui(aux2, aux2, 1);

	sem_post(&sem_den);
	return NULL;
}

/*Primeira parcela de a[k+1]*/
void *func1(void *arg){

	mpf_add_ui(aux1, y, 1);/*(1+y(k+1))*/
	mpf_pow_ui(aux1, aux1, 4);/*(1+y(k+1))^4*/
	mpf_mul(a, aux1, a);/*ak(1+y(k+1))^4*/
	
	sem_post(&sem_1);
	return NULL;
}

/*Primeira parcela de a[k+1]*/
void *func2(void *arg){
	mpf_t aux4;

	mpf_pow_ui(aux2, y, 2); /*y(k+1)^2*/
	mpf_add(aux2, aux2, y); /*y(k+1)+y(k+1)^2*/
	mpf_add_ui(aux2, aux2, 1);/*y(k+1)+y(k+1)^2+1*/
	mpf_mul(aux2, aux2, y);/*y(k+1)*(1+y(k+1)+y(k+1))^2*/
	mpf_init_set_ui(aux4, 2);
	mpf_pow_ui(aux4, aux4, (2*k+3));/*2^(2*k+3)*/
	mpf_mul(aux2, aux2, aux4);
	mpf_clear(aux4);
	sem_post(&sem_2);
	return NULL;
}

